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Abstract 

We study black hole MACHO binary formation through three-body inter- 
actions in the early universe at t ~ 10 _5 s . The probability distribution 
functions of the eccentricity and the semimajor axis of binaries as well as 
of the coalescence time are obtained assuming that the black holes are ran- 
domly formed in space. We confirm that the previous order-of-magnitude 
estimate for the binary parameters is valid within ~ 50% error. We find that 
the coalescence rate of the black hole MACHO binaries is ~ 5 x lO" 2 x 2 ±x 
events/year /galaxy taking into consideration several possible factors which 
may affect this estimate. This suggests that the event rate of coalescing bi- 
nary black holes will be at least several events per year within 15 Mpc. The 
first LIGO/VIRGO interferometers in 2001 will be able to verify whether the 
MACHOs are black holes or not. 

PACS numbers: 04.30.-w; 95.35.+d; 97.60.Lf; 98.80.-k 
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I. INTRODUCTION 



The analysis of the first 2.1 years of photometry of 8.5 million stars in the Large Mag- 
ellanic Cloud (LMC) by the MACHO collaboration suggests that the fraction 0.62±g;| 
of the halo consists of MACHOs (MAssive Compact Halo Objects) of mass O.51qJM in 
the standard spherical flat rotation halo model. The preliminary analysis of four years of 
data suggests the existence of at least four additional microlensing events with t^ur ~ 90 
days in the direction of the LMC [0. At present, we do not know what MACHOs are. 
This is especially because there are strong degeneracies in any microlensing measurements; 
the mass, velocity and distance of a lens object. The inferred mass is just the mass of red 
dwarfs. However, a tight constraint is obtained from the observations The red dwarfs 

contribute at most a few percent to the mass of the halo. 

The brown dwarfs are also restricted by the Hubble Space Telescope search. Extrapolating 
the mass function, the contribution of the brown dwarfs to the mass of the halo is less than 
a few percent ||. However, the possibility of brown dwarfs cannot be rejected if the mass 
function has a peak at the brown dwarfs since they are so dim. The possibility that the 
MACHOs are neutron stars is ruled out by the observational constraints on the metal and 
helium abundance ||. 

As for white dwarfs, assuming the Salpeter initial mass function (IMF) with an upper 
and a lower mass cutoff, the mass fraction of the white dwarfs in the halo should be less than 
10 % from the number count of the high z galaxies [0]. The observation of the chemical yield 
does not favor that the MACHOs are white dwarfs ||||. However if the IMF has a peak 
around ~ 2M , MACHOs can be white dwarfs [|TcHT2]|. Future observations of high- velocity 
white dwarfs in our solar neighborhood will make clear whether white dwarf MACHOs exist 
or not. Of course, it is still possible that an overdense clump of MACHOs exist toward the 
LMC 0. 

If the number of the high- velocity white dwarfs turns out to be large enough to explain the 
MACHOs, then star formation theory should explain why the IMF has a peak at ~ 2M Q . If it 
is not, we must consider other possibilities such that MACHOs are primordial black holes or 
boson stars. If MACHOs are black holes, however, it seems difficult to verify observationally 
whether the MACHOs are black holes or not. In fact, electromagnetic radiation from gas 
accreting to a black hole MACHO (BHMACHO) is too dim to be observed unless the velocity 
of BHMACHO is exceptionally small 0. 

Recently, however, Nakamura et al. [15] proposed the detectability of the gravitational 
wave from coalescence of BHMACHO binaries which are formed through the three body 
interaction in the early universe at t ~ 10 _5 s. The event rate of coalescing BHMACHO 
binaries was estimated as ~ 5 x 10 -2 events/year /galaxy, which suggests that we can expect 
several events per year within 15 Mpc. If this estimate is true, not only can we confirm 
whether the MACHOs are black holes or not, but also we have plenty of sources of grav- 
itational waves. However in Ref. |15| they made only order-of-magnitude arguments and 
there are uncertainties in the estimate of the event rate. Especially, the estimate of the 
semimajor and the semiminor axis of the binary formed through the three body interaction 
is not based on accurate numerical calculations. In this paper we investigate up to what 



extent the order-of-magnitude arguments in Ref. |15| are valid by calculating numerically 
the three body problem in the expanding universe. 



2 



In §2 we will review the formation scenario of the BHMACHO binaries in Ref. ||15|| . In 
§3 we will derive the basic equations of a multi-particle system in the expanding universe. 
In §4 we will show the results of the numerical calculations of the binary formation through 
the three body interaction in the expanding universe. In §5 we will obtain the probability 
distribution functions of the eccentricity and the semimajor axis of binaries as well as the 
coalescence time assuming that the black holes are randomly formed in space. We will 
estimate the event rate of coalescing binaries using this probability distribution function. In 
§6 we will consider several possible factors which may affect the estimate of the event rate, 
and we will conclude that the event rate is not reduced considerably. In §7 we will consider 
how BHMACHO binaries evolve after the binary formation. §8 will be devoted to summary 
and discussions. 

We use units of G = c = 1 in this paper. 



II. REVIEW OF BHMACHO BINARIES FORMATION SCENARIO 

We briefly review the BHMACHO binaries formation scenario in Ref. [13] to introduce 



our notations. 

For simplicity, we assume that black holes dominate the dark matter, i.e. Q = Qbhm- 
The extension to the case where black holes do not dominate the dark matter is not difficult. 
Further we assume that all black holes have the same mass and we describe it as Mbh- 
(The extension to the unequal mass case is straightforward and is given in Appendix A.) 
Primordial black holes are formed when the horizon scale is equal to the Schwarzschild radius 
of a black hole. There are some theories about the formation mechanism of primordial black 



holes p6|-|T8| . At present, however, we can not say definitely whether the black holes will 
form or not in the early universe. Ultimately, only by observational technique the existence 
of a population of primordial black holes may be established. It is therefore important to 
establish the observational signatures of primordial black holes. Our standpoint is that we 
are quite ignorant of whether large numbers of primordial black holes will form or not in the 
early universe, and it is very important to confirm the existence of primordial black holes 
observationally. Whether the results of the observation confirm the extistence of primordial 
black holes or not, it will add very much to our understanding of the universe. 
The scale factor at the time of the formation is given by 



Rf = y/M BH /H~i = 1.1 x 1(T 8 J 2 (nh*) , (2.1) 



where H eq with H' 1 = ^/Mp eq = 1.2 x 10 21 {Qh 2 ) ' cm is the Hubble parameter at the 
time of matter-radiation equality. We normalize the scale factor such that R — 1 at the 
time of matter-radiation equality. 

The mean separation of black holes x with mass Mbh at the time of matter-radiation 
equality is given by 



x = (Mbh/p, 



eq) 



,1/3 



1.2 x lO 16 (M BH /M ) 1/3 (ft/i 2 )- 4/3 cm. (2.2) 
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As a foundation for computing the distribution function of BHMACHO binaries with respect 
to the binary parameters, in Ref. [|15|] it was assumed that (i) the BHMACHOs are created 
with a distribution of comoving separations x that is uniform over the range from an initial 
physical separation equal to the black hole size to a maximum separation x = x and that 
(ii) the BHMACHOs' initial peculiar velocity is negligible compared to the Hubble flow. 
Obviously, it is more realistic to assume that the BHMACHOs are formed randomly rather 
than uniformly. We will consider the effect of the initial peculiar velocity in Section VI. D. 

Consider a pair of black holes with the same mass Mbh and a comoving separation 
x < x. These holes' masses produce a mean energy density over a sphere with the radius of 
the size of their separation as Pbh = p eq x 3 / (x 3 R 3 ) . Pbh becomes larger than the radiation 
energy density p r = p eq /R A if 

R>R m = (|) 3 . (2.3) 

After R = R m the binary decouples from the cosmic expansion and becomes a bound 
system. The tidal force from neighboring black holes gives the binary sufficiently large 
angular momentum to keep the holes from colliding with each other unless x is exceptionally 
small. 

The semimajor axis a will be proportional to xR m . Hence, we have 

4 

X 



axR m = a—, (2.4) 



x 



where a is a constant of order 0(1). To estimate the tidal torque, we assume that the tidal 
force is dominated by the black hole nearest to the binary. We denote by y the comoving 
separation of the nearest neighboring black hole from the center of mass of the binary. Then, 
from dimensional analysis, the semiminor axis b will be proportional to (tidal force) x (free 
fall time) 2 and is given by 

b ^ M BHX R (xR m r ^ p fxY^ 



(yR m ) 3 M BH 

where (3 is a constant of order 0(1). Hence, the binary's eccentricity e is given by 



\ 




(2.6) 



In Ref. [15 1, a — (3 — 1 is assumed. However, a and (3 will be different from unity so 
that calculations of the distribution functions based on an accurate estimate of a and (3 are 
necessary . This is the prime subject of the present paper. 

III. MULTI-PARTICLE SYSTEM IN THE EXPANDING UNIVERSE 

Although our main interest is in the three body problem, we formulate the problem as 
generally as possible. 
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A. Basic equations 



We treat the motion of a black hole as that of a test particle within the Newtonian 
approximation |19^pT|. We first assume that the line element is given by 

ds 2 = -(1 + 2<p)dt 2 + (1 - 2</>)i?(t) 2 dx 2 , (3.1) 

where <fi is the Newtonian potential determined by 

iT 2 A0 = Att Pbh , (3.2) 

with pBH being the energy density of black holes. For the multi-particle system, the potential 

is readily solved as </>(x) = — — — where Xj is the position of the j-th black hole. 

j i j i 
The action of the particle is given by 

Jds = J y/l + 20 - R 2 ± 2 dt ~ J (l - ^R 2 ± 2 + 0) dt. (3.3) 

Then the equation of motion is derived as 

(R 2 ±)- = -V0. (3.4) 

For the potential, 0(x;) = — — r , of the multi-particle system, the above equation 

^.Rlxi-Xjl 

is expressed as 

w ' = i^ ^ (3 - 5) 

We introduce Zi = Xj/a; and use the scale factor R as an independent variable. Then for the 
equal-mass black hole case, Eq.( [3.5|) can be written as 
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„ 1 , = M BH _ (z; - zj) 

1 i2 1 x 3 H 2 R 5 fr? Izi-z;! 

j^i i i j 



3 E ^_^l ; (3 . 6) 



87T_R fj^ |Z; — Zji 

j^l I 1 Jl 

where a prime denotes the derivative with respect to R and we have used Eq.( |2.2j ) in the 
last equality. Note that Eq. (|3.6|) does not depend on M BH . Moreover there is a scaling law, 
i.e. Eq.( |3.6| ) is invariant under the transformation defined by 

z -> Az, R -> A 3 i?, (3.7) 

where A is a constant. 
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B. Validity of Newtonian approximation 



The cosmological Newtonian approximation is valid if (1) \(f)\ <C 1 and (2) the scale 
of inhomogeneity I satisfies the relation I <C H^ 1 p0|j21| . Since the minimum separation 
between the binary black holes is a(l — e), the condition (1) is satisfied if 

M B H<o(l-e). (3.8) 

Then in terms of the initial comoving separation we have 

y/x < 5.8 x 10(x/x) 2/3 (M BH /M G )- 1/9 (nh 2 )- 2/ \ (3.9) 

where we used Eq.([2T4|) , Eq.(j2.6|) for a — (5 = 1 and the relation of (1 — e) ~ (1 — e 2 )/2. 
The condition (2) is written as 

Rx < H^R 2 . (3.10) 

Then we have 

R > 1.0 x \[)- h (x/x){M B H/M Q f^{nh 2 f^. (3.11) 
Therefore we have to choose the initial scale factor for the numerical calculations so that 



the condition Q3.ll ) is satisfied. 



IV. THREE BODY PROBLEM AND FORMATION OF BINARY BLACK HOLES 

We solve Eq.( ^.6| ) numerically for three body systems using the fifth-order Runge-Kutta 
method with the adaptive step size control |22[ . Considering the conditions for the validity of 



the Newtonian approximation derived in the previous section, we set the initial conditions 
as x/x = r/(0.1 < 7] < 1) and x = at R = lO" 3 r/(M B// /O.5M ) 1 / 3 so that Eq. (pl|) is 



satisfied, where we place a pair of black holes along the x-axis. Note that, using the scaling 
law of Eq.( |3.7| ), we see that these initial conditions are the same as x/x = Ary and x = at 
R = lO- 3 (Ar/)(M Bi? /O.5A- 6 M ) 1 /3. Therefore we can obtain the results for different M BH 
from a single numerical result. We then numerically estimate a and (3 for x and y in the 
range 0.1 < x/x < 1 and 2 < y/x < 7. The total number of the parameters we examined is 
100 for each direction of the third body. In this section we show the main results in relation 
to a and (3 first. We will discuss the dependence of a and (3 on the initial direction of the 
third body in Section VI. A. in more detail. We will also show in Section VI. D. that the 
dependence of the results on the initial conditions is small. 

In Fig.l, the trajectories of the second body (the thick curve) and the third body (the 
dotted curve) relative to the first body are shown for (a) x/x = 0.3, y/x = 2.0 and (b) 
x/x = 0.3, y/x = 4.0. 6 is chosen as 7r/4, where 6 denotes the angle between the x-direction 
and the direction of the third body. The coordinate is normalized by x. We see that the 
binary is formed through the three body interaction while the third body goes away. To see 
the accuracy of numerical calculations, we checked the time reversal of the problem. That is, 
we have re-started the numerical integration from the final time backward to the initial time. 
We have found that the differences from the true values of coordinates and velocities are 
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very small: the relative error in the coordinate position, |z ini t — z(time reversed)^ |/|z init |, 
is less than 10~ 7 and the "velocity" component, z-(time reversed), deviates from zero by at 
most 10~ 5 . 

Fig. 2 shows the semimajor axis a as a function of initial separation x/x. The filled 
triangles are numerical results. The solid line is the approximate equation, a/x = (x/x) . 
We performed the least square fitting of the numerical results assuming that a/x = a(x/x) n . 
It is found that a is well fitted by the following function 

- ~ 0.41 [-) 4.1 

x \xj 

irrespective of the direction of the third body as far as we have examined (9 = 7r/6, 7r/4, vr/3). 
The power index n is in good agreement with the analytical estimate in Eq. (|2.4| ) so that we 
will not discuss the small deviation of n from 4 from now on. 

Fig. 3 shows b/a as a function of x/y for 9 equal to (a) ir/6, (b) 7r/4 and (c) n/3. The filled 
triangles are numerical results. The solid line is the approximate equation b/a = (x/y) 3 . 
We see that the numerical results are parallel to the approximate estimate in the previous 
paper ||15|| . We performed the least square fitting of the numerical results assuming that 
b/a = (3(x/y) n . The results are given as 

3.2 

(9 = tt/6) 

3.1 

(9 = n/4) (4.2) 
(9 = vr/3) 

We see ^-dependence of (3 , which will be discussed in Section VI. A. However, as for the 
power index n, it is almost constant so that we will not discuss the small deviation of n from 
3 from now on. 

The important conclusion is that we have verified that the power dependence is in good 
agreement with the previous analytic order-of-magnitude estimate of Eq. (|2.4j ) and Eq. ( |2.5| ) 
and that numerical coefficients a and f3 are actually of order unity. In the next section we 
will assume that a and (3 are constants. For simplicity we will adopt a = 0.4 and (3 = 0.8. 




3.1 



V. PROBABILITY DISTRIBUTION FUNCTION OF BINARIES 

A. Distribution function and binary fraction 

If we assume that black holes are distributed randomly, then the probability distribution 
function P(x, y) for the initial comoving separation of the binary x and the initial comoving 
separation of the nearest neighboring black hole from the center of mass of the binary y is 

P(x, y)dxdy = ^^e- y3/si dxdy } (5.1) 
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where 0<x<y<ooso that J Q °° dx f£° dyP(x,y) = 1. Changing the variables x and y 
in Eq. (|5.1| ) to a and e with Eq.( |2.4|) and Eq.( [2.6|) , we obtain the probability distribution 
function of the eccentricity and the semimajor axis of binaries as 



f(a, e)dade = — 







a}l 2 e 



4 (ax) 3 / 2 (1 - e 2 ) 3 / 2 



exp 







(1 _ e 2)i/ 2 y a %j 



a \ 3 / 4 



dade, 



(5.2) 



where a/1 — (3 2 < e < 1, 0<a<ooso that da J^^def(a,e) = 1. For < e < 

a/1 — /3 2 , /(a, e) = 0, i.e. no such binary is formed. 

Integrating f(a,e) with respect to e in the range \/l — /3 2 < e < 1, we obtain the 
distribution function of the semimajor axis as 



fa{a)da 



3 
4 



/ a \ 3 / 4 


/ a \ 3 / 4 " 




— exp 










a 



(5.3) 



From Eq. (|5.3| ), it is found that if a = 1, the fraction of BHMACHOs that are in binaries with 
a ~ 2 x 10 14 cm and Mbh = O.5M is ~ 4% and ~ 0.4% for fi/i 2 = 1 and 0.1, respectively. 
On the other hand, if we adopt our numerical estimate of a given in the previous section, 
a = 0.4, the fraction becomes ~ 7% and ~ 0.8% for Qh 2 = 1 and 0.1, respectively. This 
estimated fraction of ~ 10 AU size BHMACHO binaries can be compared with the observed 



rate of binary MACHO events (one binary event in eight observed MACHOs) [^3|] although 
the small number statistics prevents us from stating something definite . 



B. Gravitational Waves from Coalescing BHMACHO Binaries 



We consider here short period BHMACHO binaries. Their coalescence time due to the 
emission of gravitational waves is approximately given by |24j 



* = *„(-) (i-e)" 

^ClO / 



(5.4) 



a = 2.0 x 10 



n 



' M 



BH 



Mr 



cm 



(5.5) 



where to = 10 10 year and a Q is the semimajor axis of a binary with circular orbit which 
coalesces in t - Note that Eq.( |5.4| ) is an approximation for e ~ 1 in Ref. pi ]. However it 
is also a good approximation even for e ~ 0. Eq. (|5.4| ) can be written in terms of x and y 
using Eq. (Op and Eq. (Ol) as 



X 



X 



-21 



ax\' 
a J 



7 [ — ) t . 



(5.6) 
(5.7) 



Integrating Eq. ( |5.1| ) for a given t with the aid of Eq. (|5.6|) , we obtain the probability distri- 
bution function of the coalescence time ft{t)- We should take the range of the integration 



S 



as < £ < x, x < y < oo. The first condition x < x is necessary for the binary formation. 
The second condition turns out to be [t/t) l / m x < y < (t/i)~ l l 21 x for a given t. Performing 
the integration, we have 



ft(t)dt 



3_ / A 3 / 37 
37 [tj 

3 m 3 / 37 




A 3/16 




-l/7\ 



dt 

T 



37 \t 



1 37 J t ' 



where T(x, a) is the incomplete gamma function defined by 



a) 



s x - l e 



s ds. 



(5.8) 



(5.9) 



The second equality is valid when we consider t ~ to because t^/t <^il for typical values of 
parameters, that is , t /i ~2x KT 23 for Vth 2 = 0.1, M BH = O.5M , a = 0.4 and (3 = 0.8. 

If the halo of our galaxy consists of BHMACHOs of mass ~ 0.5M Q , about 10 12 BHMA- 
CHOs exist out to the LMC The number of coalescing binary BHMACHOs with t ~ t then 
becomes ~lx 10 9 for Qh 2 = 0.1, a = 0.4 and (3 = 0.8 so that the event rate of coalescing 
binaries becomes ~ 1 x 10 _1 events/year/galaxy. This rate is slightly larger than the esti- 
mate of Ref. [15|. On the other hand, if the BHMACHOs extend up to half way to M31, the 
number of coalescing binary BHMACHOs with t ~ to can be ~ 6 x 10 9 and the event rate 
becomes ~ 6 x 10 _1 events/year/galaxy. Both of these estimates are much larger than the 
best estimate of the event rate of coalescing neutron stars based on the statistics of binary 
pulsar searches in our Galaxy, ~ 1 x 10~ 5 events/year /galaxy [25-27]. Because the first 
LIGO/VIRGO interferometers in 2001 should be able to detect BHMACHO coalescence out 
to about 15 Mpc distance, i.e., out to the Virgo Cluster ||15|| , the event rate will be several 
events per year even if we pessimistically estimate it (~ 1/100 events/year /galaxy in each 
galaxy like our own). 

In deriving the probability distribution function for the coalescence time in Eq.( |5.8| ), we 
have neglected various effects, such as the angle dependence of /3, 3-body collision, the effect 
of the fourth body, the effect of the mean fluctuation field, the initial condition dependence 
and the radiation drag. We will consider these effects in the next section. 



C. The region checked numerically 

Since we have solved three body problem only for a restricted parameter range of x and 
y, one may wonder whether our computations may not be complete. Thus we need to show 
that the parameter range of our calculations is sufficiently large. 

We have verified Eq. (|2.4j) and Eq. (|2.5| ) numerically for x and y in the range, O.lx < 
x < x, which means 

KT 3 < (f) 3 <l (5-10) 
and 2 < y/x < 7, which corresponds to 
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x\ 3 fy\ 3 oao f x 



< - < 343 - . (5.11) 

xj \xj \xj 

On the other hand, y is expressed by x and t from Eq.( |5.6|) . Therefore if we are interested 
in the coalescing binaries with the coalescence time t\ < t < t 2 , the range of y is expressed 
by 




(5.12) 

This range of y determines the probability distribution function f t (t) in Eq. fl5.8p for t\ < 
t < t 2 . In Fig. 4 the horizontal axis and the vertical axis are (x/x) 3 and exp(— (y/x) 3 ), 
respectively. The dashed lines show x = O.lx and y/x = i (i=2, 3, 4, 5, 6, 7), respectively. 
Solid lines show t x = 0.1t and t 2 = 10t for Vth 2 = 0.1, M BH = O.5M , a = 0.4 and 
(3 = 0.8 in Eq. (|5.12|) . Since the area in Fig.4 is directly proportional to the probability 
P(x,y) = d((x/x) 3 )d(e~( y / xS) ) from Eq. (|5.1| ), almost all the region we are interested in 
(O.lto ~ t < 10to) is checked numerically. So the probability distribution function f t (t) in 
Eq.flSTBD is valid for 0.1t ~ t ^ 10t though in deriving f t (t) we used Eq. (p^4|) and Eq. (p~5]) 
for the region of x and y beyond our calculations. 



VI. CONSIDERATIONS OF VARIOUS EFFECTS 

We shall consider several possible factors which may affect the estimate of the event rate 
of coalescence. 



A. Angle dependence 

So far, we have treated f3 as a constant. In reality, however, f3 has an angle dependence. 
In this subsection we investigate whether the angle dependence of (3 affects the estimate of 
the event rate. In the analytical estimate of Eq. (|2.5|) , b is proportional to the tidal force. 
Since the tidal force is proportional to sin(26 ) ), we expect that (3 is also proportional to 
sin(2#). In Fig. 5 we show the result of numerical calculations for the angle dependence of 
(3 averaged for various values of x and y/x for the exponent n = 3. We find that the angle 
dependence of (3 can be approximated by 

(3 ~ 0.8sin(2#). (6.1) 

In the previous section we have used the maximum value of j3 so that if we take this angle 
dependence into account the effective (3 will decrease. The probability distribution function 
f t {t) is proportional to /5" 21/37 since f t (t) oc t~ 3/37 oc /T 21 / 37 in Eq.(p^|) and Eq.([T7|). Hence 
qualitatively the effect of the angle dependence is to increase the event rate. 

If we consider the initial direction of the third body, the distribution function P(x, y, 6) 
for x, y and 9 is given by 

P(x,y,8)dxdyd6 = P(x,y)dxdysm6d6 = — — — e~ v ^ x dxdysm6d8, (6-2) 
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where < 9 < n/2 and we assumed that P(x, y) does not depend on 6. Integrating Eq. ( |BT2| 
for a given t with the aid of Eq. ( |5.6|) , we obtain 



ft n9 (t)dt 



'5 arcsin(i/i)? 




3/37 



3/16N 



-r 




-l/7> 



sin 



/? 7 sin 7 (2#) ( —) t = tsin 7 (2#), 
V a J 



(6.3) 
(6.4) 



where (3 is replaced with /3 sin(26 l ). The lower and upper limit of the 9 integral is determined 
by t/t{9) = 1 where the integrand in Eq.( |6.3|) becomes 0, i.e. no binary which survives up 
to t ~ t is produced. We can integrate Eq. ( |6.3j ) numerically with respect to 6 for a given 
t. For example, if we take t = t , M B h = O.5M , Vth 2 = 0.1, a = 0.4 and (3 = 0.8, then 
f t ang (t ) = 1.8 x 10- 3 /t while f t (t ) = 1.0 x 10" 3 /to- For 0.1t ^ t < 10t , we can show 
that ft ng (t ) > ft(t ). Therefore the event rate of coalescing binaries may be doubled if we 
take into account the angle dependence of (3. 



B. 3-body collision 



In deriving Eq. ([5.8|) , we take the range of the integration as < x < x, x < y < oo. 
Here we consider the range of y carefully. If y < x, the third body may be bound by the 
binary in the radiation dominated era. If the bound third body collides with the binary, 
a complicated 3-body interaction occurs. It is a difficult problem to estimate how many 
binaries whose coalescence time is ~ to are left after the complicated 3-body interaction. 
So we shall exclude such a case and estimate the minimum event rate. Namely, we shall 
restrict the range of the integration to0<x<x, x < y < oo Q. This range turns out to be 
x < y < (t/i)~ l l 2l x for a given t. Integrating Eq. (|5.1|) for a given t with the aid of Eq. (|5.6| ), 
we have 



f! bodv (t)dt 



3 

37 
3 

37 



fx 3/37 

D 

f\ 3/37 

9 



r 



r 



58 
37' 

37' 




dt 

T 



dt 

T" 



(6.5) 



The second equality is valid when we consider t ~ t . Its ratio to f t (t) in Eq. 



is 



1 The factor in front of x in the condition x < y may be more than unity for the third body 
not to be bound in this case, since the binary is more massive than a single black hole. However 
we set it unity from now on, since the factor is of order 0(1) and the conclusion that the event 
rate of coalescing binaries is not reduced so much is not changed. This statement also applies to 
Eq.(|6l2|). 
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(6.6) 



That is, the ratio of the binary formation probability without the third body being bound 
by the binary to the total binary formation probability is about 60% for t ~ to. Hence the 
3-body collision reduces the event rate at most 40%. 

If we consider the fourth body, the bound third body may not collide directly with 
the binary due to the tidal force of the fourth body to the third body. The third body 
may only act as a satellite of the binary. If this fact is taken into account, the minimum 
probability distribution function f? body (t) in Eq. (|6.5| ) increases. The semimajor axis a' and 
the semiminor axis b' of the orbit of the third body will be determined by the initial comoving 
separation of the fourth body z as 

a = a'^r (6.7) 

x 6 



and 



b' = p'(y^\' (6.8) 

respectively. In deriving these values, we treat the binary as a point mass and assume that 
the analytical estimates of Eq.( [2.4| ) and Eq. Q2.5|) are valid for this system, a' and f3' will be 
different from a and (3 respectively in this case, since the mass ratio of the total mass of the 
binary to the third body's mass is not unity. However we set a' = a and f3' = f3 from now 
on for simplicity since we are making order estimate of the effect. 

The 3-body collision will not occur when 7a < r' min is satisfied, where 7 is a constant 
which takes into account the uncertainty of the criterion for the 3-body collision. r' min is 
the minimum separation between the third body and the center of mass of the binary and 
is given by 

r'min = o'(l - e') = a' - v^^P > ^. (6.9) 

The third inequality is almost an equality because the case with a' ^> b' is considered. The 
distribution function of four bodies is given by 

27 3 -3 

Px,y,z(x, y , z) dx dy dz = —^c~ z x 2 dxy 2 dy z 2 dz. (6.10) 

To calculate the probability that the third body does not collide with the binary but is bound 
to it, the above distribution function Eq. Q6.10Q should be integrated with the constraints 

x < y < x, (6-11) 
x < z, (6.12) 

5 

z 3 < A, (6.13) 
x z 



where 5 = J(3 2 /2 / y and the last inequality comes from 7a < r' min . The condition flS.12| ) is 



necessary in order that the fourth body is not bound in the radiation dominated era. First, 
we integrate Eq. (|6.10|) with respect to z as 
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where the range of the z integral is determined by Eq. fl6.12|) and Eq. (|6.13|) . The first term 
in Eq. (j6.14| ) can be integrated for a given t with the aid of Eq. ( |5.6| ) as 



e -i d [ * \ d [ " 



[(x/S) 37 ^/*)- 21 ^ 



X 3 



X 3 



1 5-111/143(^6/143 \X 3 / \X 3 




07 / /M 348/529P 

58e I U, 



*\ 3/37 

?) 

dt 



t 



. (6. 



The upper limit of this integral is determined by Eq.( |6.11[) , and the lower limit is determined 
by Eq. ( |5.6|) and x 3 < 5y 5 /x 2 with Eq.( 6.12j ) and Eq. ( 16.13 ). The second term in Eq.( |6.14 ) 
can be integrated in the same way. 

On the other hand, the probability distribution function of the binary formation without 
the third body being bound by the binary is already obtained in Eq. (|6.5|) . Then, by summing 
the case with the third body being bound by the binary and the case without the third 
body being bound by the binary, the probability distribution function without the collision 
between the third body and the binary is obtained as 



ft body {t)dt = fi body {t)dt + A Q) 3/3? [ ^ (1 - r mA43 (0 



348/5291 \ 



m .-174/143 A 348/5291 \ v ( HI 1 

143 \t) I U43' 



V 143' Vt, 



dt 

T 



f! body (t)dt + f t {t) 



37 



58er ( f 



dt 



(6.16) 



The relative error in the second equality is about a few % for t <~ to and 7 ~ 1. The 
second term is the probability distribution function for the third body to be bound by the 
binary but not to collide with the binary. For simplicity, we set 7 = 1, which corresponds 
to 5 = 0.57 for (3 = 0.8. The ratio f? body (t )/ft(t ) is 82% for M BH = O.5M , tth 2 = 0.1, 
a = 0.4 and (3 = 0.8. Comparing this ratio to the ratio in Eq. (|6.6|) , the hierarchical three 
body bound system may be produced by about 22% of the binaries that coalesce at t ~ to, 
during the radiation dominated era. During the radiation dominated era, the probability of 
the 3-body collision is about 18% of the binaries that coalesce at t ~ to 

Because the hierarchical structure consisting of a binary and a satellite becomes unstable 
by t ~ to if r 'min/ a i s ver Y close to 1, we may have to take 7 larger than unity. Moreover 
7 will depend on the eccentricity of the binary, the inclination of the orbital plane of the 
third body and so on. Here 'unstable' means that the third body crosses the binary and 
the complicated three body interaction occurs. There are some criterions for the stability of 
the binary (e.g. Ref. PS[|). If we need to estimate the event rate accurately, we will have to 
pay attention to the value of 7. However we do not need such an accuracy here, so we set 7 
as a constant. For 7 = 10, the ratio f? body (t )/f t (t ) is 71% for M BH = O.5M , Vth 2 = 0.1, 
a = 0.4 and f3 = 0.8. 
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In addition to the stability of the hierarchical system, there may be something to be 
considered about the value of 7. If the third body is not separated enough from the binary, 
the tidal force from the third body deforms the orbit of the binary more effectively than the 
gravitational wave does. Then the estimate of the life time of the binary using Eq. fl5.4|) is 
not adequate for t ~ to, which results in the change of the event rate estimate. Note that a 
little change in the eccentricity, e, causes a large change in the life time in Eq. (|5.4|) when the 
orbit is very eccentric, 1 — e 2 1. The effect of the orbital deformation will be discussed in 
Section VII. 



C. Effect of mean fluctuation field 



In this subsection, we will estimate the tidal force from bodies other than the third body 
and how it affects the estimate of the event rate. The tidal field from the i-th BH will be 
given by 



Ti oc 



1 



x: 



3 ■ 



(6.17) 



where Xi is the comoving separation of i-th BHMACHO from the center of mass of the 



Xi is 



' 3 \ / 3 \ / 3 ' 

Xa \ , I X 5 \ I Xj 



X 3 



binary. The distribution function for x 4 , x 5 

P(x 4 , x 5 , ■ ■ ■ , x^ = e'^'^d d (J 3 

Then the mean value of (1/x 3 ) 2 for given x 3 is estimated as 

1 v <!) d(§) ■ ■ ■ Jl Js3 d($)e~^l 



(-if) 



l /-oo poo rco rco fa 



roo roo 

dV4 I 

X" Jf Jv4 



dvr, ■ ■ ■ 



dvi^i 



IVi-2 Jvi-1 Vj 



where 



x 



Thus the mean value of the tidal field for given £3 will be estimated as 

00 1 



i=3 X i 



N 

4 



foo roc roo r 

1 + f J2 / dv ± / dv 5--- dv^x \ 

j_4 Jf Jva Jvi-2 Jv 



M dVi 
2~ 



Vi-i Vf 



(6.18) 



(6.19) 



(6.20) 



(6.21) 



The first term is the contribution from the third object and the second term is that from 
the other objects. 
For i > 5, 
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(6.22) 



Thus 



i=5 



dx 1 



1 A I I!/ - 1 

/*oo roo | 

/ dyy' / dz— 

J\ Jy Z{Z — \ / 



-f* 



eS [-k dz 4^T)' 

e / Ei(-/). 



-f» 



+ / dz 



z-1 



-f* 



(6.23) 



On the other hand, for i — 4 



J 4 := e y 



00 (if 4 

/ ^4 



-t)4 



00 (if 
— < 



/-i_ e / E i(-/). 



(6.24) 



Hence 



°° /-oo /-co /-oo /•oo rim- 

j = 4 •// ^4 •''Ui-2 Jvi-i V i 



h + I = f 



-1 



(6.25) 



Therefore (tidal force by the fourth, fifth. . . objects) / (tidal force by the third object) is 
estimated as ~ /. Note that this averaged value of the tidal force by the fourth, fifth. . . 
objects can be evaluated more easily. Because the distribution is assumed to be random 
in space, the one particle distribution function for the fourth, fifth. . . objects is given by 

an uniform distribution in x > X3 with the averaged density of the BHMACHOs. Hence 
00 1 3 r°° 1 1 

y)(— ) = — zt7 / — ^4:7rx 2 dx = , . Of course, this gives the same result as before. 
~^ x\ Anx 6 Jx 3 x b x^x 6 

If y > x, i.e. / > 1, the tidal force by the fourth, fifth.... objects dominates the one by 

the third body. In deriving f t (t) in Eq.fl5.8p or ft body {t) in Eq. flBTB]) , we should take the effect 

of the fluctuation field into consideration. If the tidal force increases, (3 increases effectively, 

eventually f t (t) and ft body (t) decrease because f t (t) and ft body (t) oc t~ 3 / 37 oc (3~ 21 / 37 . To 

estimate how the minimum event rate decreases, we use 



xl 



1+1^ 

X 



.7 } 



,r 



(6.26) 
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as the tidal force for y > x. Assuming that the analytical estimate, b ~ (tidal force) x 



(free fall time) 2 , is valid as in Eq.(2J3), we have 



With this equation and Eq. ( |2.4|) , assuming that the form of a in Eq. 
Eq.fl5.4p in terms of x and y as 



. t - / x \ 37 / y 



x 



X 



(6.27) 

is valid, we write 
(6.28) 



We take the range of the integration asO<x<x, x < y < oo which means that we exclude 
the case that the third body is bound to the binary. Integrating Eq.( |5.1|) for a given t with 
the aid of Eq. (|6.28 ) in the range of x < y < 2 1 / 3 (£/T)~ 2 / 21 x, we have 



ft fluc (t)dt 



3 /A 3 / 37 



37 \i) 
3_ /£\ 3 / 37 
37 \i 



2 -21/74 



74' 



dt 

T 



V74' J t 



(6.29) 



The ratio of //' uc {t) in Eq.(|6"29D to f t (t) in Eq.Q is 



f t fluc (t) 2- 21 / 74 r (95/74,1) 



/*(*) 



T (58/37) 



0.40. 



(6.30) 



Hence the tidal force from bodies other than the third body reduces the event rate at most 
60% Q. 

D. Initial condition dependence 



So far we have assumed that the initial peculiar velocity is vanishing and the initial scale 
factor when the bodies begin to interact is R = 10 _3 (x/x)(M B ^/0.5M o ) 1 / 3 . We consider 
here whether the results crucially depend on the initial conditions or not . 

First, we consider the initial angular momentum. We have assumed that the angular 
momentum of the binary is only from the tidal force so that it is given by 



J 



2a 



lap 2 



Ml 



x 



10 



BH y^ 



(6.31) 



2 In ff UC (t), the case where the third body is bound to the binary is excluded. So we will have 
to compare it with ff body (t) not with ft(t) to evaluate only the effect of the mean fluctuation field. 
ft flUC (t)/ff b ° dv (t) = 0.67. 
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where we used Eq.( [2.4|) and Eg. ( |2.5|) . On the other hand, if the BHMACHO binary has a 
relative velocity Vf at the formation epoch, the initial angular momentum will be evaluated 
as 

J f = M BH R f xv f , (6.32) 

where Rf is given by Eq.(2.1). We note that Vf < 1. Otherwise, the BHMACHO mass 
becomes comparable with the radiation energy within the volume that the BHMACHO 
sweeps in one Hubble expansion time at the formation epoch. In this case the drag effect 
due to the radiation field will be significant so that BHMACHO will be decelerated eventually 
and Vf < 1 after all . We can now evaluate the ratio of the two angular momenta as 



J 



J 



9 x 10 




(6.33) 

where we used Eq.Q, Eq.(gJ) and Eq.fl5~f) for M BH = O.5M , Qh 2 = 0.1, a = 0.4 and 
(3 = 0.8. For t ~ t , \ Jf/J\ can be comparable to 7 for the extreme case of Vf ~ 1 since 
t /t ~ 2 x 10~ 23 . Because f t (t) oc /T 21 / 37 , the event rate becomes 30% of Eq.(|~J) if the 
initial peculiar velocity is the maximum possible value. Hence we see that the event rate is 
not reduced so much even if the initial peculiar velocity is extremely large. 

Second, we consider whether the initial peculiar velocity considerably changes a and 
the formation epoch of the binary. Consider the case that the initial peculiar velocity is 
comparable to unity when R = Rf ~ 10 -9 in Eq. fl2.1|) for Mbh ~ 0.5M Q and Qh 2 ~ 
0.1. From Eq.( |3.5| ), the peculiar velocity is damped to ~ 10~ 6 (x/a;) _1 by the time when 
R = 10~ 3 (a;/a;), since the interaction between bodies can be neglected during this period. 
Therefore it is sufficient to investigate the two body problem assuming that the peculiar 
velocity v\ is smaller than 10 _6 (x/x) _1 at R = 10~ 3 (x/x). The result is shown in Fig.6. 
Fig. 6 shows that a and the formation epoch of the binary do not crucially depend on 
the initial peculiar velocity. When Vi = — 10 _6 (x/x)~ 1 , the event rate changes because 
f t (t) oc t -3 / 37 oc a -12 / 37 . Since a is about 10% larger, the event rate is only 4% smaller. 
Note that the ratio Rx/Rx ~ 10~ 3 f j(x/x) _1 does not depend on R. 

Finally, we consider whether the initial scale factor changes a and the formation epoch 
of the binary. We calculated the two body problem using various initial scale factors. The 
results are shown in Fig. 7. As we can see from Fig. 7, a and the formation epoch of the binary 
do not strongly depend on the initial scale factor. Qualitatively the event rate increases if 
the interaction begins earlier than R = 10 _3 x/x, because a decreases. 

E. Radiation drag 

In this subsection we consider whether the force received from the background radiation 
is greater than the gravitational force between BHMACHOs or not. There are two kinds of 
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force received from the background radiation; (i) the force from the radiation through which 
a BHMACHO is traveling [^{J and (ii) the force from the radiation which a BHMACHO 



deflects, namely the dynamical friction fl30|,|3T 



First, we estimate the force from the radiation which the BHMACHO sweeps. The force 
is estimated as 

Frad ~ (radiation momentum density) x (cross section) x (velocity of the BH) 

~ gf x M 2 BH x v. (6.34) 

When R ~ R m , the ratio of the force from the radiation to the gravitational force between 
BHMACHOs is 




< 10- 13 ( - ) (6.35) 



where we set Mbh = 0.5M Q and Qh 2 = 0.1. Therefore the force from the radiation which 
the BHMACHO sweeps can be neglected for x/x ^ 10~ 3 . Since our numerical calculations 
are performed for x/x > 0.1, the radiation drag force is negligible. 

Next, we investigate whether the dynamical friction can be larger than the gravitational 
force between the BHMACHOs. When a photon passes by a BHMACHO at a distance b, 
the photon is deflected by an angle 9$ ~ AMbh/^. This deflection changes the momentum 
direction of the photon. The momentum of the photon in the incoming direction changes 
from p to p(l — cos 8d)- Thus the momentum of the BHMACHO must be changed. This 
momentum exchange causes the dynamical friction. The force that the BHMACHO receives 
due to the dynamical friction can be estimated as 

Fdyn = J (radiation momentum density) x (velocity of the BH) x (1 — cos9d)2nbdb. (6.36) 

This expression may not be precise relativistically, but this will be a good approximation 
when ?; < 1. Assuming 9d is small, 

Fdyn * / J^v9jbdb 

^^M 2 BH vlnA, (6.37) 

where In A = lnb max /b m i n is called the Coulomb logarithm. b max {b m i n ) is the maximum 
(minimum) of the impact parameter. In the cosmological situation the horizon scale sets 
the natural maximum impact parameter. In the case we are considering b m i n ~ Mbh- 
Therefore In A cannot be so large, and Fd yn is the same order as F ra d- Hence the dynamical 
friction can also be neglected. 
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VII. LIFE OF BHMACHO BINARIES AFTER FORMATION 



Finally, we consider how BHMACHO binaries evolve after the equality time. 

First, we consider the 3-body collision during the matter dominated era. For example, 
let us consider a galaxy with mass Mq ~ 1O 12 M and a radius Rq ~ lOOkpc. Then 
the density of BHMACHOs n is about n ~ l/(10 19 cm) 3 . We estimate the time scale t co u 
for a BHMACHO binary with the semimajor axis a ~ 2 x 10 15 cm (which corresponds to 
xjx ~ 0.4) to collide with other BHMACHO. Note that almost all binaries that contribute 
to the event rate have xjx ^ 0.4 as we can see in Fig.4. If we assume the velocity of the 



binary v to be the virial velocity v ~ ^Mq/Rg ~ lOOkm/s, t co u is given by 

tcou ~ — ~ 10 12 yr > 10 10 yr, (7.1) 
ncrt> 

where cr ~ 7ra 2 ~ 10 31 cm 2 is the cross section for the 3-body collision. Hence the 3-body 
collision during the matter dominated era may be small. 

Next, we consider whether the tidal force from the other BHMACHOs alters the coales- 
cence time of the binary or not. If the tidal field from the other bodies deforms the orbit 
of the binary more effectively than the gravitational wave does, the life time of the binary 
in Eq.( |5.4| ) may be different. For example, the binary can not coalesce if the increase of 
the binding energy by the tidal force is greater than the decrease by the gravitational wave 
emission. Since the binaries that contribute to the coalescence rate at present are highly 
eccentric, 1 - e 2 < 1, a little change in the eccentricity, e, causes a large change in the life 
time in Eq. (|5.4j) . If the coalescence time of the binary in Eq. (|5.4| ) is different, the probabil- 
ity distribution function for the coalescence time ft(t) in Eq.( |5.8| ) is different since we use 
Eq.( |5.4|) in deriving ft(t), and the event rate of the coalescing BHMACHO binaries may be 
reduced. 

So let us consider the tidal field from the other bodies on the orbit of the binary after 
the equality time. First, we compare the energy loss rate by the gravitational wave with the 
binding energy change rate by the tidal field from the other bodies. The average energy loss 
rate by the gravitational wave [|4[] from a binary with the eccentricity e and the semimajor 
axis a is given by 



F(GW) = 64M| H (l + je 2 + j< 
5a 5 (l-e 2 ) 7 / 2 



(7.2) 



On the other hand, the energy change rate by the tidal force from the other body is at most 
given as 

Z(tidai) _ ( tidal force) x ( velocity ) „ Mhl x ^ (7.3) 

D T B 

where Tg = a? /2Mbh is the period of the binary and D is the distance between the 
source of the tidal force and the center of mass of the binary. The ratio is 

/ D \ 3 /r\" 43 /7/\ 21 / A/for, \~ 11/6 /07) 2 ^ 22/3 
E(GW) J flitidal) I u \ I ■' \ I U \ I - u L',ll \ / 



2 x 10 25 cm/ \xj \xj \0.5M m / V 0.1 
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(7.4) 



where we used Eq. (|2.4| ) and Eq. (|2.6|) with a = 0.4, (3 = 0.8 and e ~ 1 in the first equality. 
We used Eq.( |5.6| ) in the second equality and we set t = to in the third equality. Therefore 
if D > 1 x 10 17 cm the coalescence time based on Eq. (|5.4j) is a good approximation for a 
binary with the semimajor axis a ^ 2x 10 15 cm, which corresponds to x/x < 0.4, since 
E(gw) i E(tidai) > L Thig condition is sa ti s fied in our case since the mean separation of 

3HMACH0s is ~ 2 x 10 17 cm at the equality time. 

Second, we compare the angular momentum loss rate by the gravitational wave with 
that by the tidal field of the other body. The average angular momentum loss rate by the 
gravitational wave for a binary with eccentricity e and semimajor axis a is given by 



j(GW) 



32V2M^(l + |e 2 ) 
5a 7 / 2 (l -e 2 ) 2 



(7.5) 



On the other hand, the average rate of the angular momentum change by the tidal field is 
at most given as 



j(tidal) 

Their ratio is calculated as 

j(GW) l j (tidal) 



M 2 a 2 

(tidal force) x (length) ~ g[ . 



(7.6) 
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,2\ 22/3 



xj \0.5M Q 

12/21 /^x-90/7 fy^ 
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l BH 



BH 




(7.7) 



Therefore if D > 2 x 10 20 cm 



0.4/ V°- 5M 0/ V - 1 

j{GW) j j{tidai) > i f or a binary with the semimajor axis 

o ^ 2x 10 15 cm, which corresponds to x/x ^ 0.4. If the scale factor R becomes larger 
than ~ 10 3 , the mean separation becomes larger than 2 x 10 20 , so the tidal force can be 
neglected. In the matter dominated era R is given by (t/t eq ) 2/3 , so t ~ 10 9/2 t eg at R ~ 10 3 



where t eq = ^j3/8iip eq = ^3x 3 /8ttM B h is the equality time. The increase of the angular 
momentum A J during t eq < t < 10 9 ^ 2 t eq by the tidal force can be estimated as 



IAJI 



Jdt 



l0 9 / 2 t e „ M 2 n 2 

q m BH a i. 



^ /2 ^ M 2 BH a 2 



x 3 (t/t 



-dt 



M 2 BH a 2 t. 



eq 



(7.8) 



The ratio of the increase of the angular momentum during t eq < t < 10 9 ^ 2 t eq to the initial 
angular momentum of the binary at the formation is given by 
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58 / 7 / \ 5 /21 / m2 \ 16/21 



(7.9) 



where we used Eq.(|2]4]) and Eq.( [2.5|) in the second equality. We used Eg. ([5T6|) in the third 
equality and we set a = 0.4, (3 = 0.8 and t ~ to in the last equality. Therefore if x/x < 0.4, 
which corresponds to the most binaries that coalescence at t ~ to, as can be seen from Fig.4, 
the increase of the angular momentum during t eg < t < 10 9 ^ 2 t eq can be neglected. 

To conclude, the tidal force from the other bodies can be neglected if the binary is sepa- 
rated from the other bodies by greater than the mean separation. More detailed calculations 
taking N-body effects into account are needed to confirm the above arguments on the effect 
of the tidal force of the other body on the evolution of the binary parameters after forma- 
tion. The signs of the change rate of the energy and the angular momentum in Eq. ( |7.3|) 
and Eq.( |7.6| ) are not certain so that we argued only sufficient conditions. Moreover if the 
binding energy of the binary does not change secularly but periodically under the influence 
of the tidal field, Eq. ( [7.31 ) may be an overestimate. So the effect of the tidal field may be 
weaker. 



VIII. SUMMARY 

In this paper we have discussed black hole binary formation through three body in- 
teractions in the expanding universe. We have confirmed that the order-of-magnitude ar- 
gument in Ref. |15| is valid up within an error of ~50%. Several effects have been con- 
sidered. The effect of the 3-body collision and the mean fluctuation field may reduce the 
event rate of the coalescing BHMACHO binaries about half. On the contrary the an- 
gle dependence of the tidal force may increase the event rate about twice. The results 
do not crucially depend on the initial peculiar velocity of BHMACHOs and the initial 
scale factor when the BHMACHOs begin to interact. The radiation drag does not af- 
fect the motion of BHMACHOs. After all, the probability distribution function for the 
coalescence time f t (t) in Eq.( |5.8[ ) is a good estimate. The error in the event rate esti- 
mate can be obtained by considering the minimum event rate. The minimum event rate 
can be estimated as [1 x 10 _1 (original estimate by ft(t))] x [40%(3body collision effect + 
mean fluctuation field effect)] x [30%(maximum initial peculiar velocity effect)] ~ 1.2 x 10 -2 
events/year /galaxy. Then the event rate will be 5 x 10 -2 x 2 events/year /galaxy including 
the uncertainty from the various effects in a plain fashion. This suggests that we can at least 
expect several events per year within 15 Mpc even when the event rate is minimum, 1 x 10 -2 
events/year /galaxy. This event rate of coalescing BHMACHO binaries is comparable to or 
greater than the upper limit of that of coalescing binary neutron stars |25| . The gravitational 
wave from such coalescence should be able to be detected by LIGO/VIRGO/TAMA/GEO 
network. 

We have simplified the real situation to the three body problem, so that N-body effects 
have not been fully taken into account. They are (1) the destruction of the formed binary by 
the 3-body collision between the binary and the infalling body after the equal time, (2) the 
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deformation of the orbit by the tidal field from the other body, and so on. Although these 
effects have been estimated in Section VII and the event rate estimate does not seem to be 
influenced by these effects, more detailed calculations taking N-body effects into account 
are needed to confirm this conclusion. It is possible to investigate the N-body effect by 
N-body numerical simulations. However the dynamical range of a is very large (10 5 cm < 
a < 10 16 cm), so we need to perform numerical simulations with large dynamic range using 
the biggest supercomputer. This is an important challenging numerical problem. 

Throughout this paper, we assumed that all BHMACHOs have the same mass, although 
we have outlined the extension to the unequal mass case in Appendix A. This is based on 
the assumption of a delta-function type density fluctuation at the formation. Even in this 
case of the delta-function type density fluctuation, there is a suggestion that in reality the 



IMF of primordial black holes may continue down to zero mass limit p2[. However, the IMF 



in this case has a steep rise proportional to ~ M 3 at the lower mass end and an exponential 
cut-off near the horizon mass. Hence the picture of the delta-function-like IMF seems to 
be valid. However, in the case of a general spectrum of the density fluctuation, we should 
consider binaries made from different mass BHMACHOs. 

Although we have assumed that the initial distribution of BHMACHOs is random, the 
high density region may have a strong correlation. Presumably this depends on the black 
hole formation process or the initial density perturbation spectrum. If a strong correlation 
existed, more binary BHMACHOs may be formed. This is also an interesting future problem. 
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APPENDIX A: EXTENSION TO UNEQUAL MASS CASE 

In this appendix, we estimate the semimajor axis and semiminor axis of the BHMACHO 
binary with unequal mass. We only give the order-of-magnitude estimate along the line of 
Ref. |L5| and Section II. 



We describe the mass function of black holes as F(M), which is normalized as 
Jq°° F(M)dM = 1. The average mass of black holes Mbh can be obtained as Mbh = 
Jq 00 MF(M)dM. The mean separation of black holes at the time of matter-radiation equal- 
ity is given by x = (M BH / p eg ) 1//3 , where we assumed that the average in space is equal to the 
ensemble average. Consider a pair of black holes with masses Mi and M 2 and a comoving 
separation x. This pair will decouple from the cosmic expansion if its mean energy density 
Pbh = (Mi + M 2 )/(2a; 3 i? 3 ) becomes larger than the radiation energy density p r = p eq jR^. 
In terms of R, this condition can be written as 




R>Rm^\^hr ~ =t - , (Ai) 
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where £ = 2Mbh/{M\ + M 2 ). The semimajor axis a will be proportional to xR m and is 
given by 



2M BH \ c „x 4 
M 1 + M 2 J x 3 ^ x 3 



« = a(Tr^Trl3r=^r 5 . ( A2 ) 



where a is a constant of order 0(1). Consider a black hole with mass M 3 in the nearest 
neighborhood of the binary. Let its comoving separation from the center of mass of the 
binary be y. Then the semiminor axis b will be proportional to (tidal force) x (free fall 
time) 2 and is given by 

6 = a ' 3 { W^f) { (mTMw ) = {wTmJ e UJ = '" 9 UJ a - (A3) 

where $ is a constant of order 0(1) and rj = 2M 3 /(M 1 + M 2 ). a and j3 may depend on mass. 

If we assume that black holes are formed randomly, then the probability distribution 
function P(x, y, M 1 ,M 2 , M 3 ) is 



9t 2 ?/ 2 1 1-1 

P(x, y, Afi, M 2 , M 3 )dxdydM 1 dM 2 dM 3 = ——^—e~ y ' dxdyF(Mi)F(M 2 )F(M 3 )dM 1 dM 2 dM 3 , 



(A4) 



where we assumed that x, y do not depend on mass. Eq.(|5.4|) can be written in terms of x 
and y using Eq .(|A2|) and Eq . (|A3|) as 



/ \ '17 / \ —21 

~ f x\ 61 I y 



t = t {* I ' (A5) 



t = (r ] (3) 7 ( — \ t . (A6) 
l«o/ 



y 

~_\ 4 



Integrating Eq. (|A4j ) for a given t with the aid of Eq. (|Aq) , we obtain the probability distri- 
bution function of the coalescence time ft neq (t) for the unequal mass case. We should take 
the range of the integration as < x < £ _1 / 3 :r, x < y < 00. The first condition x < £ _1 / 3 x 
is necessary for the binary formation so that R m < 1. The second condition turns out to be 
(t/ty/ 16 x < y < £~ 37 / 63 (t/t)~ 1//21 x for a given t. Performing the integration, we have 



poo poo poo Q / + \ 3/37 

f^m=l I I |(D 



dt 

T 



x F(Mi)F(M 2 )F(M 3 )dMidM 2 dM 3 . (A7) 

To integrate the above equation with mass, we need accurate values of a and j3, as well as 
assuming the form of the mass fuction F(M). This is left as our future problem. 



23 



REFERENCES 



C. Alcock et al, Astrophys. J. 486, 697 (1997). 
M. Pratt, in Proc. 3d Microlensing Workshop, 1997 (unpublished). 
J. N. Bahcall, C. Flynn, A. Gould, and S. Kirhakos, Astrophys. J. Lett. 435, L51 (1994). 

C. Flynn, A. Gould, and J. N. Bahcall, Astrophys. J. Lett. 466, L55 (1996). 

D. S. Graff and K. Freese, Astrophys. J. Lett. 456, L49 (1996); D. S. Graff and K. Freese, 
Astrophys. J. Lett. 467, L65 (1996). 

D. Ryu, K. A. Olive, and J. Silk, Astrophys. J. 353, 81 (1990). 
S. Chariot and J. Silk, Astrophys. J. 445, 124 (1995). 
B. K. Gibson and J. R. Mould, Astrophys. J. 482, 98 (1997). 
R. Canal, J. Isern, and P. Ruiz-Lapuente, Astrophys. J. Lett. 488, L35 (1997). 
G. Chabrier, L. Segretain, and D. Mera, Astrophys. J. Lett. 468, L21 (1996). 
F. C. Adams and G. Laughlin, Astrophys. J. 468, 586 (1996). 

B. D. Fields, G. J. Mathews, and D. N. Schramm, Astrophys. J. 483, 625 (1997). 
T. Nakamura, Y. Kan-ya, and R. Nishi, Astrophys. J. Lett. 473, L99 (1996). 
Y. Fujita, S. Inoue, T. Nakamura, T. Manmoto, and K. E. Nakamura, Astrophys. J. Lett, 
(to be published). 

T. Nakamura, M. Sasaki, T. Tanaka, and K. S. Thorne, Astrophys. J. Lett. 487, L139 
(1997). 

J. Yokoyama, Astron. Astrophys. 318, 673 (1997). 

M. Kawasaki, N. Sugiyama, and T. Yanagida, |hep-ph/9710259| (unpublished). 
K. Jedamzik, Phys. Rev. D 55, 5871 (1997). 

P. J. E. Peebles, The Large-Scale Structure of the Universe (Princeton University Press, 
Princeton, 1980). 

T. Futamase, Mon. Not. R. Astron. Soc. 237, 187; Phys. Rev. D 53, 681 (1996); T. Fu- 
tamase and M. Sasaki, Phys. Rev. D 40, 2502 (1989). 
M. Shibata and H. Asada, Prog. Theor. Phys. 94, 11 (1995). 

W. H. Press et al, Numerical Recipes in Fortran 2nd Ed. (Cambridge University Press, 
1992). 

D. P. Bennett et al, |astro-ph/9606012] (unpublished). 
P. C. Peters, Phys. Rev. B 136, 1224 (1964). 

E. S. Phinney, Astrophys. J. Lett. 380, L17 (1991). 

R. Narayan, T. Piran, and A. Shemi, Astrophys. J. Lett. 379, L17 (1991). 
E. P. J. van den Heuvel and D. R. Lorimer, Mon. Not. R. Astron. Soc. 283, L37 (1996). 
V. Szebehely and K. Zare, Astron. Astrophys. 58, 145 (1977). 

C. Hogan, Mon. Not. R. Astron. Soc. 185, 889 (1978). 
S. Chandrasekhar, Astrophys. J. 97, 255 (1943). 

J. Binney and S. Tremaine, Galactic Dynamics (Princeton University Press, Princeton, 
1980). 

[32] J. C. Niemeyer and K. Jedamzik, |astro-ph / 9709072] (unpublished). 



24 



FIGURE CAPTION 



Fig.l. The trajectories of the second body (thick curve) and the third body (dotted curve) 
relative to the first body for (a) x/x = 0.3, y/x = 2.0 and (b) x/x = 0.3, y/x = 4.0. 
9 = 7r/4. The coordinate is normalized by x. 

Fig. 2. The semimajor axis a as a function of initial separation x/x. The filled triangles are 
numerical data. The solid line is the approximate equation a/x = (x/x) 4 . 

Fig. 3. 6(semiminor axis)/a(semimajor axis) as a function of x/y for (a) 9 = tt/6, (b) 7r/4 
and (c) tt/3. The filled triangles are numerical data. The solid line is the approximate 
equation b/a = (x/y) 3 . 

Fig. 4. The region we have checked numerically, i.e. x = O.lx in Eq.( |5.10| ) and y/x = i 
(i=2,3,4,5,6,7) in Eq.( |5.11| ), and the region that corresponds to O.lio <t< 10t . The 
horizontal axis is scaled as (x/x) 3 and the vertical axis is scaled as exp(— (y/x) 3 ) so 
that the area in the figure is directly proportional to the probability. We can see 
almost all the region we are interested in is within the range that we have checked 
numerically for O.lto < t < 10t . 

Fig. 5. The angle dependence of j3 assuming that the functional form of b is as in Eq.( p.5|) 
is shown. j3 has an angle dependence as (3 oc sin(2#). 9 is the angle between the line 
that connects the binary and the line that connects the third body and the center of 
the binary. 

Fig. 6. This figure illustrates whether the initial peculiar velocity crucially changes a or the 
formation epoch of the binary. The evolution of the relative distance between the 
binary with different initial peculiar velocity at Ri = 10 _3 (x/x) is shown. The case of 
Vi = 10~ 6 (x/x) _1 is the upper one. The case of Vi = is the middle one. The case of 
Vi = — 10 _6 (x/x) _1 is the lower one. We can see the dependence is week. 

Fig. 7. This illustrates whether the initial scale factor when the bodies begin to interact 
crucially changes a or the formation epoch of the binary. The cases of Ri = 10~ 3 (x/x), 
Ri = 10 _4 (x/x) and Ri = 10~ 5 (x/x) are shown. We can see the dependence is weak. 
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